home *** CD-ROM | disk | FTP | other *** search
/ PD ROM 1 / PD ROM Volume I - Macintosh Software from BMUG (1988).iso / Programming / Programming Tools / FORTRAN Routines / PREDICT.FOR < prev    next >
Encoding:
Text File  |  1984-12-30  |  1.8 KB  |  76 lines  |  [TEXT/ttxt]
  1. $NOFLOATCALLS
  2. $NODEBUG
  3. $STORAGE:2
  4. c***********************************************************************
  5.  
  6.     subroutine rkp(x,ys,h,neqn,imeth,kt,est,yold,ynew,w,func)
  7.  
  8.  
  9.     implicit double precision (a-h,k,p-z)
  10.     external func
  11.     dimension w(13,neqn),kt(neqn),ys(neqn),est(neqn)
  12.     dimension al(12),b(12,12),yold(neqn),ynew(neqn)
  13.     common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
  14.      &               a12,a13,er1,er3,er4,er5,er6,er7,er8,
  15.      &               er9,er10,er11,er12,er13,inum
  16.     xs=x
  17.     do 5 i=1,neqn
  18.        ys(i)=yold(i)
  19. 5    continue
  20.     call func(xs,ys,kt,neqn)
  21.     do 7 i=1,neqn
  22.        w(1,i)=kt(i)
  23. 7    continue
  24.  
  25.     do 20 j=2,inum
  26.        j1=j-1
  27.        xs=x+h*al(j1)
  28.        do 15 i=1,neqn
  29.           ksum=0.d0
  30.           do 10 m=1,j1
  31.          ksum=ksum+b(m,j1)*w(m,i)
  32. 10          continue
  33.           ys(i)=yold(i)+h*ksum
  34. 15       continue
  35.        call func(xs,ys,kt,neqn)
  36.        do 17 i=1,neqn
  37.           w(j,i)=kt(i)
  38. 17       continue
  39. 20    continue
  40.  
  41. c****    evaluate YNEW[i] and the error estimates EST[i]
  42.  
  43.     if(imeth.EQ.1) then
  44.     do 30 i=1,neqn
  45.        ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
  46.      &                         +a5*w(5,i)+a6*w(6,i))
  47.        est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
  48.      &                  +er6*w(6,i))
  49. 30    continue
  50.     end if
  51.  
  52.     if(imeth.EQ.2) then
  53.     do 40 i=1,neqn
  54.        ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
  55.      &                         +a5*w(5,i)+a7*w(7,i)
  56.      &                         +a8*w(8,i))
  57.        est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
  58.      &                  +er6*w(6,i)+er7*w(7,i)+er8*w(8,i))
  59. 40    continue
  60.     end if
  61.  
  62.     if(imeth.EQ.3) then
  63.     do 50 i=1,neqn
  64.        ynew(i)=yold(i)+h*(a1*w(1,i)+a6*w(6,i)+a7*w(7,i)
  65.      &                 +a8*w(8,i)+a9*w(9,i)
  66.      &                 +a12*w(12,i)+a13*w(13,i))
  67.        est(i)=h*(er1*w(1,i)
  68.      &            +er6*w(6,i)+er7*w(7,i)+er8*w(8,i)+er9*w(9,i)
  69.      &            +er10*w(10,i)+er11*w(11,i)+er12*w(12,i)
  70.      &            +er13*w(13,i))
  71. 50    continue
  72.     end if
  73.  
  74.     return
  75.     end
  76.